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ABSTRACT 

Bright  Field  (BF)  electron  tomography  (ET)  has  been  widely  used  in  the  life  sciences  to  characterize  biological 
specimens  in  3D.  While  BF-ET  is  the  dominant  modality  in  the  life  sciences  it  has  been  generally  avoided  in  the 
physical  sciences  due  to  anomalous  measurements  in  the  data  due  to  a  phenomenon  called  “Bragg  scatter”  -  vis¬ 
ible  when  crystalline  samples  are  imaged.  These  measurements  cause  undesirable  artifacts  in  the  reconstruction 
when  the  typical  algorithms  such  as  Filtered  Back  Projection  (FBP)  and  Simultaneous  Iterative  Reconstruction 
Technique  (SIRT)  are  applied  to  the  data.  Model  based  iterative  reconstruction  (MBIR)  provides  a  powerful 
framework  for  tomographic  reconstruction  that  incorporates  a  model  for  data  acquisition,  noise  in  the  mea¬ 
surement  and  a  model  for  the  object  to  obtain  reconstructions  that  are  qualitatively  superior  and  quantitatively 
accurate.  In  this  paper  we  present  a  novel  MBIR  algorithm  for  BF-ET  which  accounts  for  the  presence  of  anoma¬ 
lous  measurements  from  Bragg  scatter  in  the  data  during  the  iterative  reconstruction.  Our  method  accounts  for 
the  anomalies  by  formulating  the  reconstruction  as  minimizing  a  cost  function  which  rejects  measurements  that 
deviate  significantly  from  the  typical  Beer’s  law  model  widely  assumed  for  BF-ET.  Results  on  simulated  as  well 
as  real  data  show  that  our  method  can  dramatically  improve  the  reconstructions  compared  to  FBP  and  MBIR 
without  anomaly  rejection,  suppressing  the  artifacts  due  to  the  Bragg  anomalies. 

1.  INTRODUCTION 

Bright  Field  (BF)  electron  tomography  (ET)  has  been  widely  used  in  the  life  sciences  to  characterize  biological 
specimens  in  in  3D.1  BF-ET  typically  involves  acquiring  microscope  images  (of  transmitted  electrons)  corre¬ 
sponding  to  various  tilts  of  a  sample,  and  using  an  algorithm  on  the  acquired  “tilt-series”  to  reconstruct  the 
attenuation  coefficient  of  the  object.  In  most  cases  due  to  the  geometry  of  the  acquisition  and  mechanical 
limitations  of  the  tilting  stages,  BF-ET  is  a  limited  angle  parallel  beam  tomography  modality. 

While  BF-ET  is  the  dominant  modality  in  the  life  sciences  it  has  been  generally  avoided  in  the  physical 
sciences  due  to  contrast  reversals2  from  Bragg  scatter  in  crystalline  samples.  Bragg  scatter  occurs  when  the 
crystal  lattice  is  oriented  in  such  a  manner  that  the  incident  electrons  are  elastically  scattered  away  from  the 
direct  path  leading  to  an  anomalous  measurement  uncharacteristic  of  attenuation  due  to  thickness  alone.  The 
presence  of  Bragg  anomalies  in  the  data  can  result  in  artifacts  since  typical  tomographic  reconstruction  algorithms 
(like  FBP  and  SIRT3)  do  not  account  for  these  effects. 

Model  based  iterative  reconstruction  (MBIR)  provides  a  powerful  framework  for  tomographic  reconstruction 
that  incorporates  a  model  for  data  acquisition,  measurement  noise  and  for  the  object  to  obtain  reconstructions 
that  are  qualitatively  superior  and  quantitatively  accurate. 4-6  While  Levine7  has  applied  MBIR  to  BF-ET  in 

Further  author  information: 

S.V.  Venkatakrishnan:  svenkata@purdue.edu 

C.A.  Bouman:  bouman@purdue.edu 

J.P.  Simmons:  Jeff.Simmons@wpafb.af.mil 

L. F.  Drummy:  Lawrence.Drummy@wpafb.af.mil 

M.  De  Graef:  mdg@andrew.cmu.edu 


Computational  Imaging  XI,  edited  by  Charles  A.  Bouman,  Ilya  Poliak,  Patrick  J.  Wolfe, 
Proc.  of  SPIE-IS&T  Electronic  Imaging,  SPIE  Vol.  8657,  86570A  -©2013  SPIE-IS&T 
CCC  code:  0277-786X/1 3/$18  -  doi:  10.111 7/1 2.201 3228 


SPIE-IS&T/ Vol.  8657  86570A-1 

Distribution  A.  Approved  for  public  release;  distribution  unlimited. 

Downloaded  From:  http://spiedigitallibrary.org/  on  11/12/2014  Terms  of  Use:  http://spiedl.org/terms 


Tilt: -10°  Tilt :  0°  Tilt:  +10° 


Bragg  scatter 


Figure  1.  Illustration  of  the  effect  of  Bragg  scatter  on  a  real  TEM  data  set  of  Aluminum  nanoparticles.  The  figure  shows 
BF  images  corresponding  to  three  different  tilts  of  the  specimen.  Note  that  certain  spheres  turn  dark  (fewer  counts)  and 
then  again  turn  bright.  Due  to  the  orientation  of  the  crystals,  electrons  are  scattered  away  from  the  BF  detector  leading 
to  fewer  electrons  being  collected. 

the  case  of  thick  specimens  and  shown  that  it  can  improve  the  quality  of  the  reconstructions,  his  work  deals  with 
cases  where  there  is  no  anomalous  Bragg  scatter  in  the  measurement. 

In  this  paper,  we  present  a  MBIR  algorithm  for  BF-ET  which  identifies  and  rejects  the  anomalous  mea¬ 
surements  in  the  data  during  the  reconstruction.  We  use  a  Beer’s  law  forward  model  and  combine  it  with  a 
prior  model  for  the  material  to  formulate  the  reconstruction  as  minimizing  a  cost  function.  The  cost  function 
is  designed  so  that  it  rejects  the  measurements  that  deviate  significantly  from  the  assumed  forward  model  due 
to  Bragg  scatter.  We  then  develop  a  fast  algorithm  to  minimize  the  cost  function.  Our  algorithm,  which  is 
based  on  the  iterative  coordinate  descent  (ICD),  works  by  constructing  a  substitute  to  the  original  cost4  at  every 
point,  and  minimizing  this  new  function.  This  operation  significantly  lowers  the  computational  complexity  of 
the  optimization  and  speeds  up  the  overall  convergence  of  the  algorithm.  Intuitively,  our  method  starts  with  a 
initial  reconstruction,  forward  projects  and  compares  it  with  the  true  measurements  and  in  case  the  deviation 
is  large,  classifies  those  measurements  as  anomalous.  Using  the  new  set  of  non-anomalous  measurements,  the 
sample  is  reconstructed  and  we  repeat  the  process.  Thus  as  the  reconstruction  progresses,  measurements  are 
constantly  being  monitored  and  the  anomalous  ones  are  rejected. 

We  apply  our  method  to  simulated  data  sets  with  Bragg  scatter  like  anomalies.  Results  show  that  our 
method  can  significantly  improve  the  reconstructions  compared  to  FBP  and  MBIR  without  anomaly  rejection, 
suppressing  the  artifacts  that  arise  due  to  the  anomalous  measurements.  We  also  apply  our  method  to  a  real 
BF  data  set  and  demonstrate  that  it  can  reduce  the  artifacts  compared  to  FBP  and  MBIR  without  anomaly 
rejection. 

The  organization  of  the  rest  of  the  paper  is  as  follows.  In  section  2  we  introduce  a  statistical  model  for 
the  measurement,  combine  it  with  a  prior  model  for  the  material  and  formulate  the  MBIR  cost  function  which 
accounts  for  anomalous  Bragg  measurements.  In  section  3  we  develop  an  efficient  algorithm  to  minimize  the 
cost  function.  In  section  4  we  present  results  from  a  simulated  data  set,  followed  by  results  from  a  real  data  set. 
Finally  in  section  5  we  draw  our  conclusions. 

2.  MEASUREMENT  MODEL  AND  COST  FORMULATION 

The  goal  of  BF-ET  is  to  reconstruct  the  attenuation  coefficient  at  every  point  in  the  sample.  An  electron  beam  is 
focused  on  the  material  and  the  electrons  that  are  transmitted  through  the  sample  are  captured  by  a  BF  detector 
to  obtain  a  single  image.  The  sample  is  then  tilted  along  a  fixed  axis  and  the  process  is  repeated.  Thus  at  the 
end  of  the  acquisition  we  obtain  a  collection  of  BF  images  which  can  be  used  for  tomographic  reconstruction  of 
the  attenuation  coefficients. 
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The  reconstruction  in  the  MBIR  framework  is  typically  given  by  the  joint-MAP8  estimate 


(/,  </>)  =  argmin  {-  log p{g\f,  (/>)  -  log p(f)}  (1) 

f,4> 

where  g  is  the  vector  of  measurements,  /  is  the  vector  of  unknown  voxels  (attenuation  coefficients),  0  is  a  vector 
of  unknown  parameters,  p(g\f,  <j>)  is  the  likelihood  and  p(f)  is  the  prior  probability  for  the  unknown  voxels.  Next 
we  derive  the  above  cost  function  for  the  case  of  BF-ET  accounting  for  Bragg  scatter  in  the  measurements. 

First  we  present  an  expression  for  the  likelihood  term,  assuming  there  is  no  Bragg  scatter  in  the  measurement. 
Let  A i  be  the  electron  counts  corresponding  to  the  ith  measurement  and  Ad  be  the  counts  that  would  be  measured 
in  the  absence  of  the  sample.  We  model  the  attenuation  of  the  beam  through  the  material  using  Beer’s  law. 
Thus  the  projection  integral  corresponding  to  the  *th  measurement  is  given  by  log  There  can  be  cases 

in  which  the  dosage  Ad  is  not  measured  and  we  can  include  it  as  a  unknown  nuisance  parameter  in  the  MBIR 
framework.  If  g  is  a  M  x  1  vector  with  gi  =  —  log(Aj),  /  is  a  N  x  1  vector  of  unknown  attenuation  coefficients 
of  the  material,  d  =  —  log(AD)  is  an  unknown  offset,  then  using  a  Taylor  series  approximation  to  the  likelihood 
function,  formed  by  assuming  A,;’s  are  independent  Poisson  random  variables,9  gives 


-logp(ff|/,d)  «  ^\\g-  Af-dl\\\  +  h(g)  (2) 

1  M 

=  2  “  Ai’*f  -  rf)2  A«  +  MsO 

z  i=  1 

where  A  is  a  M  x  N  forward  projection  matrix,  A  is  a  diagonal  matrix  with  entries  An  set  to  be  inversely 
proportional  to  the  variance  of  the  measurement  and  h(.)  is  some  function  of  the  data.  Ignoring  the  effect  of 
electronic  noise,  we  set  A„;  =  Aj.4  We  note  that  our  formulation  can  account  for  more  sophisticated  models  as 
introduced  in,10  but  in  this  paper  we  focus  on  using  the  Beer’s  law  as  it  has  been  found  to  be  accurate  for  a 
certain  class  of  thin  (<  Ifixn)  samples. 

The  above  likelihood  can  be  directly  used  to  formulate  the  MBIR  cost  function  for  BF-ET  when  measurements 
are  consistent  with  our  model.  However  due  to  Bragg  scatter,  some  of  the  measurements  are  anomalous,  i.e.  they 
deviate  significantly  from  the  model.  Bragg  scatter  typically  results  in  electrons  being  scattered  away  from  the 
direct  path  resulting  in  fewer  counts  than  would  be  expected.  Fig.  1  shows  an  example  of  three  tilts  from  a  BF 
tilt  series  with  regions  having  significant  Bragg  scatter  (indicated  using  an  arrow).  A  precise  way  of  accounting 
for  this  would  be  to  model  the  mechanisms  that  causes  Bragg  scatter,  but  this  can  be  very  complicated;  and  so 
in  this  work  we  account  for  Bragg  anomalies  by  rejecting  those  measurements.  We  use  the  penalty  function 


Pt(x) 


x2  |x|  <  T 
T2  |x|  >  T 


to  limit  the  influence  of  anomalous  measurements.  This  function  (Fig.  2)  plays  a  similar  role  to  the  weak-spring 
potential11  used  to  model  image  priors,  where  it  is  used  to  limit  the  influence  of  pixels  across  an  edge.  Thus  the 
new  likelihood  is  given  by 

1  M 

-  log  p(g\f,d)  =  iS9i  “  Ai’*f  -  rf)  (v^))  +  h(g)  (3) 

i= 1 

where  At  *  is  the  *th  row  of  the  forward  projection  matrix  A,  h(g)  is  a  constant  independent  of  /  and  d.  The 
threshold  T  can  be  set  so  that  for  the  Bragg  scattered  measurements  the  ratio  of  the  data  fit  error,  (gi  —  A^f  —  d), 
to  the  noise  standard  deviation  in  the  measurement,  (-^==),  is  greater  than  T.  T  can  be  set  as  a  user  input.  Note 
that  the  above  likelihood  does  not  correspond  to  a  proper  density  function  since  the  area  under  the  corresponding 
density  function  is  unbounded.  However  this  formulation  can  still  be  used  to  compute  the  joint  MAP  estimate. 
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Figure  2.  Illustration  of  the  penalty  function  /3t  used  for  the  likelihood  term  with  T  =  2.  Large  model  mismatch  errors 
are  penalized  by  restricting  their  influence  on  the  overall  cost  function. 


To  model  the  prior  p(f)  we  use  a  q-generalized  Gaussian  Markov  random  field  (qGGMRF)12  of  the  form 


P(f) 

s(f) 


P(fj  ~  fk) 


tf“  exp  {— s(/)} 

H  wjkp(fj  ~  fk) 
{i,fe}e  M 


f j  fk 

°7 

Q 

c  + 

f j  fk 

°7 

q-p 

(4) 


where  Zf  is  a  normalizing  constant,  M  is  the  set  of  pairs  of  neighboring  voxels  (e.g.  a  26  point  neighborhood), 
and  p,  q,  c  and  07  are  qGGMRF  parameters.  The  weights  Wjk  are  inversely  proportional  to  the  distance  between 
voxels  j  and  k,  normalized  to  1.  Typically  l<p<g<2is  used  to  ensure  convexity  of  the  function  p(.),  thereby 
simplifying  the  subsequent  MAP  optimization.  We  fix  q  =  2  and  c  =  0.001  so  that  the  prior  behaves  similar  to  a 
GGMRF,13  and  it  has  a  bounded  second  derivative  which  is  a  useful  property  for  the  subsequent  optimization. 

Substituting  (3)  and  (4)  into  (1),  the  reconstruction  is  obtained  by  minimizing  the  cost 


c(/,  d) 


1 

2  (( 9i  -  Ai,*f  -  d)^/K^j  +  8(f). 


Alternately  we  can  define 


({Ji  d)  Aj2  I  (pi 

T2  \(g  i 


A,*f  ~  d)y/Au\  <  T 
Aq*/  -  d)v/A)i'|  >  T 


(5) 


so  the  cost  function  can  be  written  as 


c(f>d)  =  ^2,PT,i(f,d)  +s(f)  (6) 

Additionally  we  will  constrain  /  >  0  as  it  is  physically  meaningful  to  have  positive  values  of  the  attenuation 
coefficients.  Thus  the  MBIR  BF-ET  reconstruction  is  given  by 


(id) 


argmin  c(/,  d) 
f>o,d 


SPIE-IS&T/Vol.  8657  86570A-4 

Distribution  A.  Approved  for  public  release;  distribution  unlimited. 

Downloaded  From:  http://spiedigitallibrary.org/  on  11/12/2014  Terms  of  Use:  http://spiedl.org/terms 


3.  OPTIMIZATION  ALGORITHM 


The  cost  function  (6)  is  non-convex,  and  thus  finding  the  global  minimum  is  difficult.  Here  we  attempt  to  find  a 
desirable  local  minimum  of  the  cost.  We  minimize  the  cost  function  in  (6)  using  the  ICD14  algorithm.  In  ICD  we 
start  with  a  initial  value  for  the  variables,  and  then  they  are  updated  one  at  a  time;  so  that  each  update  lowers 
the  value  of  the  cost  function.  Minimizing  (6)  with  respect  to  each  variable  can  be  computationally  expensive 
due  to  complicated  form  of  the  likelihood  and  prior  terms,  so  we  can  instead  construct  a  substitute  to  the  original 
function  and  minimize  this  new  function.  A  substitute  function  to  the  original  function  is  constructed  so  that  it 
bounds  the  original  function  from  above,  and  so  that  minimizing  the  substitute  function  results  in  a  lower  value 
of  the  original  cost. 

Our  goal  is  to  find  a  substitute  function  to  the  original  cost  so  that  it  can  be  easily  minimized  with  respect 
to  each  voxel  or  the  unknown  parameter  d.  We  find  a  substitute  function  for  the  likelihood  terms  and  the  prior 
terms  separately  and  then  sum  them  to  form  a  single  substitute  function  for  the  original  cost. 

We  design  substitute  functions  Qr,i(f  ,d;  f  ,d')  for  each  of  /3 T,»(/>d)  in  (6)  at  a  given  point  (/',d').  In 
particular 


Qr,i{f,  d;  f' ,d!) 


(, 9i  -  A,*/  ~  d)2Au 

rp2 


-d')VKi\<T 

l(Si-A,./'-d'y^|  >t 


is  a  substitute  function  for  /3r,i(-)  as  shown  in  Appendix  A. 

We  can  find  a  substitute  function  for  each  potential  function  p(fj  —  fk)  at  the  point  f  of  the  form 


(7) 


Pifi-fk-jj-fk)  =  +  (8) 

Using  such  a  form  results  in  a  simple  closed  form  updates  for  the  voxels  during  the  optimization.  The  values  of 
a,jk  and  bjk  can  be  derived  as  shown  in  Appendix  B  and  are  given  by 


Ujk 

(  P  (fj  fk)  f!  _L  f! 

=  I  (fj-fi) 

(9) 

K(o)  f'j  =  f'k 

bjk 

(10) 

Thus  a  substitute  function  to  s(f)  at  / 

=  f  is 

«(/;/')  = 

wjkP(fj  -  fk;  f'j  -  fk )• 

(ll) 

{jM  eW 


3.1  Algorithm 

Based  on  the  present  value  of  (f,d)  which  we  denote  ( f,d ')  we  define  the  following  indicator  variable: 


bi  = 


I  (Si 
I  (ff* 


A,*/' 

A,*/' 


d')VKi\  <  t 
d')VKi\  >  T 


(12) 


Intuitively  bi  indicates  which  measurements  are  classified  as  anomalous  and  which  are  not,  based  on  the  current 
state  of  the  reconstruction.  Using  (7)  and  (11)  a  substitute  function  to  the  original  cost  (6)  at  (f  ,d')  is 

1  M 

QtU,  d;  f,  d')  =  -Y,  QtM ,  d;  f,  d')  +  *(/;  /') 

2—1 

1  M  M 

=  -J2(9i-Ai,J-d)2Aiibi+  M/;-/Ar/D+2E(1-yr2  (13) 

i=i  U,fc}eAA  i—i 
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function  [f,d\  4-  Reconstruct^,  /',  d',  R) 

%Inputs:  Measurements  g ,  Initial  reconstruction  /',  Initial  dosage  d! ,  Fraction  of  entries  to  reject  R 
%Outputs:  Reconstruction  /  and  dosage  parameter  d 
e'  =  g  —  Af  —  d!  1 

r  =  0  >  Initial  fraction  to  reject 

T  ■£-  oo,  b  <—  1  >  Initially  do  not  reject  any  measurement, 

while  Stopping  criteria  is  not  met  do 

for  each  voxel  j  in  random  order  do  o  Iterate  over  all  voxels. 

M 

@2  4—  Aabi 

i= 1 
M 

6\  i - e(A nbi(Ajj) 


i—  1 

for  k  £  A fj  do 

Compute  substitute  function  parameter  ajk  using  (9) 

end  for 

wjkajkfk+^2f'j-Si 

u*  < - ^ - 

}  j  Wjkajk+S 2 
fceA/j 

fj  4—  max(u* ,  0) 
e'  4-  e'  -  ( fj  -  f'j)A*j 

fj  fj  _ 

Update  b  using  (12)  >  Computing  the  new  Qt  function 

end  for 

d!  4—  Update  d  using  (16)  >  Dosage  parameter  update 

Update  e' 

Update  b  using  (12) 

If  r  <  R,  then  r  =  r  +  R/ 10  >  Increment  the  rejection  threshold 

Compute  new  T  >  Sort  the  array  of  e*  *  V An  and  set  T  =  rth  percentile 

Update  b  using  (12) 

end  while 

/  /)  d  ■*—  d 

end  function 


Figure  3.  MBIR  algorithm  for  BF  data  with  Bragg  scatter.  The  algorithm  works  by  constructing  a  substitute  to  the 
original  function  based  on  the  current  values  of  the  voxels  and  dosage  parameter  and  minimizing  this  substitute  function 
with  respect  to  a  single  variable.  The  process  is  then  repeated.  The  algorithm  can  be  efficiently  implemented  by  keeping 
track  of  the  error  sinogram.15  Further  the  rejection  ratio  r  is  progressively  increased  till  it  reaches  the  target  value  R  to 
prevent  the  algorithm  from  getting  stuck  in  undesirable  local  minima. 
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3.1.1  Voxel  Update 

The  voxels  are  updated  one  at  a  time  in  random  order  similar  to4  in  order  to  speed  up  the  overall  convergence 
of  the  algorithm.  To  minimize  with  respect  to  voxel  j .  we  can  fix  ff  =  f'k  \/k  €  {1,  •  •  •  ,  M}  \  {j}  and  d  =  d!  in 
(13).  The  cost  function  we  need  to  minimize  is 


Csub (u)  =  dm  +  y  (u  -  /j)2  +  Y  wjkP(u  ~  f'k,  fj  ~  f'k)- 

keNj 

where  An  =  Anbj,  §i  =  — (e')1  A (A*  j),  d2  =  (A*  j)*A(^4*  j),  j  is  the  jth  column  of  the  forward  projection 
matrix  A  ,  e'  =  g  —  Af  —  d'  1  ,  /j  is  the  present  value  of  voxel  j,  and  A fj  is  the  set  of  all  neighbors  of  voxel  j. 
Since  p(u  —  /(;  /'  —  f'k)  is  quadratic  in  u,  the  minimum  of  csub(w)  has  a  closed  form  and  is  given  by 


Y  wjkajkf'k  +  6*2 fj  -  <5*1 

fee  A/j 

^  '  WjkCtjk  T  d‘2 
feeA/j 


(14) 


Enforcing  the  positivity  constraint,  the  update  for  the  voxel  is 

fj  4—  max  (u* ,  0) 


(15) 


3.1.2  Nuisance  Parameter  Update 

In  order  to  minimize  the  substitute  function  with  respect  to  the  dosage  parameter  d, 
the  substitute  function  (13)  Qt{/' ,d\  f  ,d')  with  respect  to  d  and  set  it  to  zero.  This 
for  d  as 

(e'YAl 

a  4 —  — -—z — 

1*  A1 

d  4—  d  +  a  (16) 


we  take  the  derivative  of 
gives  the  optimal  update 


3.1.3  Multi-resolution  Initialization 

The  optimization  can  be  further  sped  up  using  a  multi-resolution  initialization. 16  In  multi-resolution  initializa¬ 
tion,  we  perform  a  reconstruction  at  a  coarser  resolution  (larger  voxel  sizes)  and  use  its  output  to  initialize  a 
finer  resolution  reconstruction.  This  transfers  the  computational  load  to  the  coarser  scale  where  the  optimization 
can  be  done  quickly  due  to  to  the  reduced  dimensionality  of  the  problem.  Since  our  prior  behaves  similar  to  a 
GGMRF,13  we  adapt  the  scaling  parameter  af  according  to  Eq.28  in.17 

We  set  the  rejection  threshold  parameter  T  in  the  algorithm  indirectly,  via  an  user  input  R,  the  approximate 
fraction  of  the  total  measurements  that  the  are  affected  by  Bragg  scatter.  We  refer  to  this  as  the  target  rejection 
rate.  Given  a  value  of  R ,  T  can  be  chosen  so  that  RM  measurements  are  not  used  in  the  cost  function. 

The  algorithm  is  terminated  if  the  ratio  of  the  average  change  in  the  magnitude  of  the  reconstruction  to 
the  average  magnitude  of  the  reconstruction  is  less  than  a  preset  threshold.  Further,  in  order  to  prevent  the 
algorithm  from  getting  stuck  in  undesirable  local  minima,  the  rejection  percentage  is  gradually  increased  to  the 
desired  target  rejection  rate  (and  therefore  T  is  gradually  decreased).  The  MBIR  BF-ET  algorithm  for  a  single 
resolution  is  summarized  in  Fig.  3. 


4.  RESULTS 


In  this  section  we  compare  three  algorithms  for  BF-ET  -  FBP,  MBIR  without  Bragg  anomaly  correction  and 
MBIR  with  anomaly  correction.  While  BF-ET  has  been  avoided  in  the  physical  sciences  because  algorithms  like 
FBP  are  know  to  be  unreliable  for  this  type  of  data,  we  include  it  to  indicate  how  an  algorithm  designed  for  the 
application  (MBIR)  can  outperform  a  standard  algorithm  (FBP).  FBP  is  performed  in  Matlab  using  the  iradon 
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Figure  4.  Simulated  BF  data  corresponding  to  a  phantom  of  spheres  for  three  successive  tilts.  The  figure  in  the  center 
shows  the  simulated  Bragg  scatter  obtained  by  lowering  the  counts  by  about  50%.  This  is  an  anomaly  in  the  projection 
data,  and  if  not  accounted  for,  can  cause  artifacts  in  the  reconstruction. 


x  —  z  slice  from  phantom 


Figure  5.  A  single  x  —  z  slice  from  the  phantom  data  used  for  qualitative  comparison  against  various  reconstruction 
algorithms. 


(a)  FBP  (b)  MBIR  without  Bragg  correction  (c)  MBIR  with  Bragg  correction 

Figure  6.  Comparison  of  BF  reconstructions  on  a  data  set  with  few  tilts  having  anomalous  measurements,  (a)  shows 
a  single  x  —  z  slice  from  a  FBP  reconstruction,  (b)  shows  the  MBIR  reconstruction  without  Bragg  correction.  The 
reconstruction  is  comparable  to  the  phantom  because  the  fraction  of  anomalous  measurements  is  relatively  low.  The 
strong  streaking  artifacts  are  significantly  suppressed  compared  to  FBP.  (c)  shows  the  reconstruction  with  Bragg  anomaly 
correction  with  the  rejection  threshold  set  to  5%  .  The  reconstruction  is  superior  to  the  case  in  which  we  apply  no 
correction  as  well  as  FBP.  All  images  are  scaled  in  the  range  of  0  —  7.45  x  10”3  nm_1. 
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(a)  FBP  (b)  MBIR  without  Bragg  correction  (c)  MBIR  with  Bragg  correction 

Figure  7.  Comparison  of  BF  reconstructions  for  a  data  set  with  a  high  percentage  of  the  tilts  having  Bragg  scatter,  (a) 
shows  a  single  x  —  z  slice  from  a  FBP  reconstruction,  (b)  shows  the  MBIR  reconstruction  without  Bragg  correction.  The 
reconstruction  has  streaks  because  of  the  Bragg  scatter  but  much  lesser  compared  to  FBP.  (c)  shows  the  reconstruction 
with  Bragg  anomaly  correction  with  the  rejection  threshold  set  to  10%.  The  method  effectively  suppresses  the  artifacts 
in  (a)  and  (b),  and  produces  a  more  accurate  reconstruction.  All  images  are  scaled  in  the  range  of  0  —  7.45  x  10~3  nm-1. 


command.  The  output  is  clipped  to  be  positive.  For  the  MBIR  reconstructions,  <jf  and  R  are  chosen  to  obtain 
the  best  visual  quality  of  reconstruction.  The  offset  d  is  initialized  by  taking  the  mean  value  of  the  signal  in  a 
void  region  of  the  sample. 

For  the  first  experiment  we  use  a  3-D  phantom  containing  spheres  with  an  attenuation  coefficient  of  7.45  x  10-3 
nm-1.  The  sample  has  a  dimensions  of  256  nm  x  512  nm  x  512  nm  (z  —  x  —  y  respectively).  The  phantom  is 
forward  projected  at  141  tilts  in  the  range  of  —70°  to  +70°  in  steps  of  1°  with  a  dosage  A u  =  1850  counts  using 
the  Beer’s  law  model  with  Poisson  noise.  At  certain  tilts  («  40%)  some  of  the  measurements  are  decreased  by  50% 
to  simulate  Bragg  scatter  like  effects  (Fig.  4).  While  accurate  simulation  of  Bragg  scatter  is  very  complicated,  we 
attempt  to  demonstrate  some  artifacts  that  can  occur  even  in  this  very  simple  case  and  show  how  our  algorithm 
can  handle  it.  We  reconstruct  a  3D  volume  consisting  of  12  (x  —  z )  slices. 

Fig.  5  shows  a  single  x  —  z  slice  from  the  original  phantom.  Fig.  6  shows  the  reconstructions  of  the  corre¬ 
sponding  x  —  z  slice.  The  FBP  reconstruction  (Fig.  6  (a))  has  streaking  artifacts  due  to  Bragg  scatter  as  well 
as  the  absence  of  a  prior  model  for  the  object.  The  MBIR  reconstructions  (Fig.  6  (b)  and  (c))  significantly 
suppress  the  streaking  artifacts  compared  to  FBP.  We  observe  that  the  MBIR  result  without  Bragg  rejection 
(Fig.  6(b))  is  qualitatively  comparable  to  the  MBIR  with  the  Bragg  rejection  (Fig.  6(c))  though  it  has  some 
streaking  artifacts.  This  suggests  that  when  the  fraction  of  anomalous  measurements  is  low,  the  algorithms  with 
and  without  the  Bragg  rejection  produce  qualitatively  comparable  results.  However  Table  1  shows  that  MBIR 
with  the  Bragg  anomaly  correction  produces  a  quantitatively  more  accurate  reconstruction. 

Fig.  7  shows  the  x  —  z  slice  reconstructed  using  the  different  algorithms  when  we  use  only  a  subset  of  47 
tilts  from  the  phantom  data  set.  Most  of  these  tilts  have  Bragg  anomalies  and  hence  the  fraction  of  anomalous 
measurements  is  much  higher  in  this  data  set.  The  FBP  reconstruction  (Fig.  7(a))  has  strong  streaking  artifacts 
indicating  why  it  has  not  been  used  for  BF  reconstructions  in  the  physical  sciences.  The  MBIR  without  Bragg 
anomaly  correction  (Fig.  7  (b))  shows  prominent  streaking  artifacts  in  the  reconstruction  even  though  it  is 
much  lesser  than  in  FBP.  However  MBIR  with  anomaly  correction  (Fig.  7(c))  produces  a  reconstruction  which 
effectively  suppresses  these  artifacts.  Table  1  shows  that  MBIR  with  the  Bragg  anomaly  correction  significantly 
improves  the  quantitative  accuracy  of  the  reconstruction  compared  to  FBP  as  well  as  MBIR  without  anomaly 
correction. 

Fig.  8  shows  a  x  —  z  slice  («  581  nm  x  900  nm)  reconstructed  from  a  real  sample  of  spherical  Aluminum 
nanoparticles.  The  BF-TEM  data  consists  of  15  tilts  in  the  range  of  —70°  to  +70°  in  steps  of  10°.  The  FBP 
reconstruction  (Fig.  8  (a)),  has  strong  streaking  artifacts.  The  reconstruction  using  the  MBIR  algorithm  with  no 
anomaly  correction  (Fig.  8  (b)),  also  has  streaking  artifacts  similar  to  those  in  the  simulated  data  set  of  Fig.  7 
but  much  lesser  than  in  the  case  of  FBP.  Fig.  8  (c)  shows  that  using  the  Bragg  anomaly  correction  can  result  in 
reconstructions  which  have  significantly  reduced  streaking  artifacts  compared  to  MBIR  without  Bragg  anomaly 
rejection  (Fig.  8  (b))  as  well  as  FBP  (Fig.  8  (a)).  The  Bragg  rejection  threshold  was  set  to  reject  10%  of  the 
data. 
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Table  1.  Comparison  of  the  Root  Mean  Square  Error  of  the  reconstruction  with  respect  to  the  original  phantom  for  various 
scenarios.  MBIR  with  Bragg  anomaly  correction  produces  quantitatively  more  accurate  reconstructions. 


Data  Set 

Algorithm 

Bragg  Correction 

R.MSE  (nm-1) 

Limited  Bragg  (Fig.  6) 

FBP 

No 

13  x  1CT4 

Limited  Bragg  (Fig.  6) 

MBIR 

No 

5.7  x  1CT4 

Limited  Bragg  (Fig.  6) 

MBIR 

Yes 

4.87  x  1CT4 

More  Bragg  (Fig.  7) 

FBP 

No 

23  x  1CT4 

More  Bragg  (Fig.  7) 

MBIR 

No 

11.95  x  10"4 

More  Bragg  (Fig.  7) 

MBIR 

Yes 

6.52  x  1CT4 

(a)  FBP  (b)  MBIR  without  Bragg  correction  (c)  MBIR  with  Bragg  correction 

Figure  8.  A  single  x  —  z  slice  reconstructed  from  a  BF-TEM  data  set  of  Aluminum  sphere  nanoparticles.  The  FBP 
reconstruction  (a)  has  very  strong  streaking  artifacts,  suggesting  why  it  has  been  avoided  for  BF-ET.  The  MBIR  algorithm 
with  the  anomaly  rejection  (c)  is  superior  to  the  case  in  which  we  apply  no  correction,  suppressing  the  streaking  artifacts 
seen  in  (b).  In  the  case  of  MBIR,  the  circular  cross  section  of  the  spherical  particles  are  clearly  visible  compared  to  FBP. 
All  images  are  scaled  in  the  range  of  0  —  4.0  x  10~3  nm-1  and  the  rejection  threshold  is  set  to  10%. 
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5.  CONCLUSION 


In  this  work  we  presented  a  MBIR  algorithm  for  BF-ET  which  can  significantly  decrease  the  artifacts  in  the 
reconstruction  due  to  anomalous  Bragg  scatter.  Our  method  works  by  modeling  the  image  formation  and 
the  sample  being  imaged  to  formulate  a  cost  function  which  rejects  measurements  that  do  not  fit  the  model 
accurately  as  a  part  of  the  reconstruction.  Results  on  simulated  and  real  data  sets  demonstrate  that  our  method 
can  effectively  suppress  the  artifacts  due  to  Bragg  scatter,  producing  qualitatively  and  quantitatively  accurate 
reconstructions. 


APPENDIX  A.  SUBSTITUTE  FUNCTION  FOR  LIKELIHOOD 


Theorem  A.l.  Each  QT,i(f,d;f',d')  is  a  substitute  function  to  f)T,i(f,d). 

Proof.  Let  us  define  the  following  functions.  QiK— l  R,  /3  :  R  — ►  R  such  that 


P(x) 


x 2  \x\  <  T 

T2  \x\  >  T 


Q(x\x’) 


x2  \x’\  <  T 
T2  \x’\  >  T 


Then  Q(x\x’)  is  a  substitute  function  to  /3(x).  If  we  define  a  function  hi  :  Rw+1  — >  R  such  that  hi(f,d)  = 
(. gt  —  Ai^f  —  d)y/Kii.  Then  Q(hi(f,d)-,hi(f',d'))  is  a  substitute  function  to  /3(hi(f,d))  by  Property  7.6  in.18 
Thus  the  theorem  is  proved  by  recognizing  Qr,i(f,d;  f  ,d')  =  Q(hi(f,d);hi(f',d'))  and  j3T,i{f,d)  =  P(hi{f,d)). 
□ 


APPENDIX  B.  SUBSTITUTE  FUNCTION  FOR  PRIOR 

In  order  to  find  a  suitable  substitute  function  to  the  prior  s(f)  each  of  the  potential  functions  p(fj  —  fk)  can  be 
replaced  by  a  function  p(fj  —  fk-  ft  —  f'k)  which  satisfy  the  following  properties4 

Pifi-fkif’j-f'k)  >  P(fj-fk)  V/jER  (17) 

p'Ui-fcfj-fk)  =  p'U'j-fk)  (18) 

where  f  is  the  point  of  approximation.  Intuitively  (17)  ensures  that  the  substitute  function  upper  bounds 
the  original  potential  function  and  (18)  ensures  that  the  derivatives  of  the  original  function  and  the  substitute 
function  are  matched  at  the  point  of  approximation.  We  use  a  substitute  function  of  the  form 

p(fj-fk-J'-K)  =  ^(fj-ftf  +  bjk  (19) 

because  it  results  in  a  simple  closed  form  update  for  a  given  voxel.  Thus  we  need  to  find  the  values  of  ajk  and 
bjk  which  satisfies  (17)  and  (18).  Taking  the  derivative  of  the  substitute  function  (19)  and  matching  it  to  the 
derivative  of  the  original  potential  function  we  get 

7^  fk 
f'j  =  fk 

To  choose  bjk  we  set  the  value  of  the  original  potential  function  and  substitute  function  to  be  the  same  at 
the  point  of  approximation  ft.  This  gives 

bjk  =  P(fj  ~  f'k)  ~  -  fk)2 


p' (f'j- f'k) 

(fi-f'k) 
P"(  0) 
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